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ABSTRACT 

Given  any  analytic  non-characteristic  initial  value 
problem  it  is  possible  to  obtain  a  truncated  power  series 
as  an  approximation  to  the  solution.   The  complete  power 
series,  from  which  the  truncation  comes,  may  not  have  a 
domain  of  convergence  large  enough  to  cover  the  area  of 
interest.   Two  methods  for  continuing  these  truncated 
power  series,  similar  to  analytic  continuation  for  a 
complete  power  series,  are  given.   Proofs  are  also  given 
that  these  two  methods  of  continuation  converge,  as  the 
number  of  terms  in  the  truncated  power  series  goes  to 
infinity,  to  the  same  values  as  the  analytic  continuation 
of  the  complete  power  series.   These  two  methods  are 
applied  to  the  problem  of  finding  the  supersonic  flow 
past  a  blunt  axisymmetric  body  and  the  results  are 
compared. 


-  2  - 


Nyo-2881 

TABLE  OP  CONTENTS 

Page 

Abstract  2 

Section 

1.  Introduction l\. 

2.  The  First  Method  of  Continuation 12 

3.  The  Second  Method  of  Continuation 23 

[(.,   Results  of  the  Test  Calculations •  3^4- 

Appendix  A.   Details  of  the  Text  Problem  .  .  [\3 
Appendix  B.   The  Arithmetic  Used  in 

the  Test  .     '. 59 

Table  1 62 

Table  2 63 

Table  3 65 

Table  1^ 6? 

References 68 


-  3  - 


TWO  METHODS  USING  POWER  SERIES  FOR  SOLVING 
ANALYTIC  INITIAL  VALUE  PROBLEMS 

1.  Introduction. 

Over  a  period  of  years  there  has  been  a  growing 
Interest  in  problems  in  partial  differential  equations 
which  are  improperly  posed  (in  the  sense  of  Hadamard). 
Much  of  classical  physics  lead  to  problems  in  differential 
equations  which  are  well  posed  such  as  the  classical 
Initial  and  boundary  value  problems.   These  well  posed 
problems  often  admit  of  an  elegant  theoretical  treatment 
as  well  as  practical  computational  methods.   This  is 
not  true  of  most  of  the  problems  which  are  improperly 
posed.   These  problems  are  often  difficult  to  treat 
mathematically.   Special  cases  have  been  treated  by 
C.  Pucci,  M.  Picone,  F.  John  and  others,  but  no  general 
methods  seem  available,   (See  [l],  [2],  [3J,  and  [I4.].) 
Many  of  these  problems  fall  in  the  category  of  analytic 
non-characteristic  initial  value  problems  which  are  not 
hyperbolic.   Even  though  the  Cauchy-Kowalewski  Theorem 
guarantees  that  there  exists  one  and  only  one  analytic 
solution  for  such  problems  (at  least,  in  the  small),  it 
does  not  guarantee  that  other  solutions  do  not  exist, 
and  even  the  analytic  solutions  need  not  depend 
continuously  on  the  initial  data.   This  paper  is  only 
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concerned,  however,  with  finding  the  analytic  solution 
to  analytic  non-characteristic  initial  value  problems; 
non-analytic  solutions,  which  probably  cannot  exist  for 
such  problems,  are  not  considered. 

One  problem  from  fluid  dynamics  which,  stated 
in  an  appropriate  way,  can  be  placed  in  this  category 
is  the  supersonic  flow  past  a  blunt  axisymmetric  body. 
As  is  well  known,  a  detached  shock  forms  in  front  of 
the  body.   Ideally,  it  would  be  desirable  to  give  the 
shape  of  the  body  and  from  this  to  determine  the  shape 
of  the  shock  and  the  various  flow  quantities  (fluid 
velocity,  pressiire,  density,  etc  )  between  it  and  the 
bodyo   However,  this  is  not  an  initial  value  problem 
and  furthermore  presents  great  mathematical  difficultieso 
On  the  other  h-nd  if  the  shape  of  an  axisymmetric  shock 
is  given  as  an  analytic  function  then  the  flow  quantities 
may  be  determined  as  analytic  functions  on  the  shock 
surface  through  the  use  of  the  Rankin e-Hugoniot  shock 
conditions o   The  equations  of  motion  coupled  with  these 
conditions  then  presents  an  analytic  initial  value 
problemo   The  general  characteristics  of  this  flow  are 
well  known  if  a  reasonable  shape  for  the  shock  is  assumed, 
The  equations  change  type    depending  on  the  velocity  of 
the  flowo   The  flow  is  subsonic  in  a  region  in  front 
of  the  tip  of  the  body  and  the  equations  are  elliptic; 
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the  flow  becomes  supersonic  away  from  the  tip  and  the 
equations  become  hyperbolic;  see  Figure  1,    page   6  „ 
The  problem  is,  therefore,  of  mixed  type,  and  is  probably 
improperly  posed  in  the  sense  that  continuous  analytic 
pertixrbations  of  the  shock  shape  will  not  produce 
continuous  changes  in  the  shape  of  the  body.   This 
may  well  be  due  to  the  concept  of  "improperly  posed. " 
Since  this  problem  comes  from  physics,  it  is  probably 
"properly  posed"  if  a  more  restricted  concept  is  used 
(eog.  only  solutions  which  have  an  a  priori  bound  in 
a  given  dom.ain  are  compared  with  the  original  solution)  o 
All  things  considered,  this  problem  makes  an  excellent 
test  case  for  any  computational  method  which  claims 
to  solve  analytic  initial  value  problems  in  general. 
The  general  characteristics  are  known,  the  constancy  of 
entropy  alonr  streamlines  and  the  constancy  of  the 
Bernoulli  number  afford  exacting  checks,  and  the  problem 
is  of  practical  as  well  as  theoretical  interesto 

The  first  method  which  suggests  itself  in  solving 
any  analytic  non-characteristic  initial  value  problem 
is  the  use  of  power  serieso   Following  the  method  used 
in  the  proof  of  the  Cauchy-Kowalewski  Theorem,  it  is 
possible  to  compute  all   the  low  order  coefficients  of 
the  solution  expressed  as  a  power  series  about  any 
point  on  the  initial  surfaceo   This  truncated  power 
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series  should  yield  good  approximations  to  the  values 
of  the  solution  at  least  neaj?  the  point  of  expansion. 
This  method  was  first  applied  to  the  detached  shock 
problem  by  Lin  and  Rubinov  (see  [5])  and  later  by 
several  others.   This  expansion  was  pursued  furthest 
by  R.  D,  Rlchtmyer  (see  [6]).   He  was  the  first  to  devise 
systematic  recursion  relations  for  the  coefficients 
adequate  for  the  use  on  a  large-scale  computer.   (These 
relations  are  described  in  Appendix  A.)   His  interest 
was  primarily  to  investigate  the  use  of  power  series 
for  such  problems  as  opposed  to  the  interest  of  the  fluld- 
dynamlcists.   He  placed  the  problem  on  the  Univac, 
carried  nine  figures  of  accuracy,  and  computed  all 
coefficients  of  the  flow  quantities  up  to  order  llj.. 
The  results  were  excellent  near  the  point  of  expansion, 
but  there  was  som.e  question  as  to  whether  the  position 
of  the  body  was  accurately  computed.   To  answer  this 
question  it  was  decided  to  recede  the  problem  for  the 
IBM  7OI4.  with  increased  accuracy.   This  report,  then, 
is  the  direct  outgrowth  of  that  decision.   In  some 
sense,  the  present  report  is  a  continuation  of 
R.  D.  Rlchtmyer 's;  but,  for  the  convenience  of  the 
reader,  it  was  decided  to  make  the  two  reports 
independent. 

The  problem  was  receded  to  carry  Sh   bits  of  accuracy. 
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and  to  handle  coefficients  up  to  order  19.   It  then 
became  clear  that  the  power  series  of  the  flow 
quantities  diverged  before  the  body  was  reached. 
This  difficulty  was  predicted  by  M.  D.  Van  Dyke 
(see  [7J)«   Therefore,  it  became  clear  that  a  method 
for  continuing  truncated  power  series  was  needed  similar 
to  analytic  continuation  for  the  complete  power  series. 
The  obvious  thing  to  do  is  to  treat  the  truncated 
power  series  as  if  it  were  the  complete  power  series 
and  to  analytically  continue  this  truncated  power  series. 
Clearly,  this  is  doomed  to  fail.   The  truncated  power 
series  is  but  a  polynomial,  say  in  x  of  order  k,  and 
analytic  continuation  about  a  point  x   of  this  polynomial 
simply  expresses  this  polynomial  in  powers  of  (x-x  ) 
and  does  not  change  its  values  anywhere.   As  k  goes  to 
infinity  (i.e.  more  and  more  terms  of  the  original 
power  series  are  taken),  the  resulting  polynomials 
must  diverge  for  any  point  which  is  neither  in  nor  on 
the  boundary  of  the  dom.ain  of  convergence  of  the  original 
power  series.   This  is  merely  a  restatement  of  the  fact 
that  a  power  series  must  diverge  for  any  point  which 
is  strictly  outside  of  its  domain  of  convergence. 
However,  if  a  certain  proportion  of  the  high  order  terras 
of  the  polynomials  expressed  in  powers  of  (x-x  )  are, 
discarded,  the  approximations  will  converge  to  the 
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desired  value  as  k  goes  to  infinltyo   This  is  called 
the  first  method  of  continuation,  is  described  in 
detail  in  Section  II,  and  a  proof  of  convergence  as  k 
goes  to  infinity  is  giveno 

The  second  method  of  continuation  described  in 
detail  in  Section  III  uses  a  different  approach.   As 
in  the  first  method,  a  truncated  power  series  is  computed 
on  the  initial  surface;  but,  now,  this  truncated  power 
series  is  used  to  approximate  initial  data  on  a  new 
initial  s-urface,,   Prom  this  new  Initial  data  a  new 
truncated  power  series  is  computed  using  the  differential 
equations  again.   If  the  new  initial  surface  satisfies 
certain  restrictions  and  if  the  order  of  the  new 
truncated  power  series  is  a  certain  fraction  of  the 
order  of  the  old  truncated  power  series,  then  it  is 
shown  in  Section  III  that  this  new  truncated  power  series 
converges  as  the  number  of  terms  goes  to  infinity  for 
arguments  outside  of  the  domain  of  the  original  power 
serieSo   Surprisingly  enough,  the  second  method  was  the 
first  discovered.   Various  methods  were  discussed,  and 
finally  the  second  method  was  tried  on  the  machine. 
The  results  were  extremely  good.   In  attempting  to 
justify  these  results  the  much  simpler  first  method  was 
hit  upon. 

Section  IV  describes  and  compares  the  results 
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when  the  two  methods  were  applied  to  the  detached 
shock  problem.   Appendix  A  describes  the  detached  shock 
problem  in  detail  and  gives  the  details  of  the  actual 
test  computations,  and  Appendix  B  describes  the 
■unnormalized  double-precision  floating-point  arithmetic 
used. 

The  author  would  like  to  express  his  gratitude  to 
Dr.  R.  D.  Richtmyer,  Prof.  E.  Isaacson,  and 
Prof.  P.  John  whose  continued  interest  and  helpful 
suggestions  made  this  paper  possible. 
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2,  The  First  Method  of  Continuation. 

Description  and  proof  of  convergence.   As  was  already 
stated,  it  is  possible  to  compute  the  low  order 
coefficients  of  the  analytic  solution  of  an  analytic 
initial  value  problem.   (The  actual  computation  of 
these  coefficients  is  discussed  in  detail  in  Appendix  A.) 
In  this  section  the  solution  will  be  thought  of  as  an 
analytic  function  of  which  the  low  order  coefficients 
of  its  power  series  are  known  and  no  fiirther  reference 
to  the  analytic  Initial  value  problem  will  be  made. 

Let  P(x. )  be  an  analytic  function  of  the  variables 
x-,,..o,x  which  is  expressible  as  a  power  series  about 
the  origin.   Let  D  be  the  domain  of  convergence  of 
this  power  series  (that  is,  the  largest  open  set  which 
contains  only  points  for  which  the  power  series  is 
absolutely  convergent).   Let  x.  be  any  fixed  point 
in  D  «   Under  the  assumptions  P(x. )  may  be  expressed 
as  a  power  series  about  the  point  x. .   Let  D,  be  the 
domain  of  convergence  of  this  new  power  series,  and 
let  X.  be  any  fixed  point  in  D,  ,   The  problem,  then, 
is  to  find  an  approximation  to  P(x. )  when  only  the 
truncated  power  series  about  the  origin  of  F(x. )  is 
known  o 

Introduce  the  new  coordinate  system  £..    =   x.-x.. 
Let  f (C^)  =  P(43_+x°)  and  let  ^°  =  x^-x°,  then  P(x]^)=f(C?) 
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Assuming  only  the  terms  of  order  less  than  or  equal  to 
k  of  P(x. )  expanded  about  the  origin  are  known,  let 
P,  (x,  )  be  this  truncated  power  series,  and  let 
f,  (4.)  =  P,  (4.+X?).   Pinally  let  f  A^,)   be  the  terras 
of  f,  (^. )  of  order  less  than  or  equal  to  [sk]  where  s 
is  some  constant  between  zero  and  one  ([skj  is  the 
greatest  integer  in  sk).   That  is,  the  truncated  power 
series  F,{x. )  is  written  in  terms  of  the  new  coordinate 
system  E,. ,    and  then  only  terms  of  £,.    of  order  less  than 
or  equal  to  [sk]  are  retained  to  obtain  the  final 
appr  oximat  ion. 

Theorem:   Under  the  above  assumptions,  there 

exists  an  s,  0  <  s  <  1,  such  that  f   ,  (£?)  ->  P(x. )  as 
J  J  s,k  ^1        i 

k  ->  00  , 

Before  the  proof  it  is  convenient  to  consider 
the  properties  of  the  function  g(z), 

g(z)  =  AB^     ^ 


Z/t   i1-z  * 
z  (l-z; 

for  0  <  s  <  1  where  A  and  B  are  positive  constants. 

The  minima  of  this  function  occur  at  the  extremes  of 

the  range:   g(z)  =  A  when  z  =  0,  g(z)  =  AB  when  z  =  1. 

The  maximum,  g(z)  =  A(B+1),  occurs  when  z  =  B/(1+B). 

Thus  the  equation,  g(z)  =  1,  has  at  most  two  roots, 

see  Figure  2,  page   15« 

Choose   R.   such  that  R.  is  in  D   and  |x?I  <  R. . 
1  1        o      '  i  '     1 
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This  is  always  possible,  see  [9].   Let  A  =  max  |x?|/R. 

£  X     2. 

and  B  =  max    |.g°|/|x°I.      Choose    an  s,    0  <   s  <  1,    such  that 


g(z) 


(0,AH> 


(B/(1+B),A(1+B|) 
I. 


(0,1)^ 


/ 


(1,AB) 


/ 


<—  s- 


(0,0) 


I         (1,0) 


Figixre    2 


s  is  smaller  than  the  roots  of  g(z)  =  1.   Under  this 
assumption  it  will  be  established  in  the  proof  that 
f  ,  (5°)  ->  F(x5")  as  k  ->  CO  .   Since  A  <  1,  it  is 

S ,  iC   X  X 

always    possible    to    choose    an   s  which   satisfies    the 
above   property. 

Proof:      Since   F(x.  )    is    assiJined  to   be    analytic   at 
the    origin,    it   has   a  power   series   expansion: 

00     n,  n 

F(x.  )   =  >       5  a^  X,  ...   x^      . 

1  -n.    — x        J nT,...,n        1  q 

n=0  nT+...+n^-n        1*         *    q  ^ 

1  q 

(All  subscripts  are  assumed  to  be  non-negative.) 

F(x. )  may  therefore  be  written  in  the  following  way: 
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n  CO  m^  ^m 

P(xJ)    =IZIZ  ^      ...   m     «1  •••    <   ^ 

m=0  in-+..«+in  =m        1'***'    q  ^ 

1  q 

where     b  equals 

1*         '    q 

1  q  ^  ^ 

(    (^)    =   0   if  n  <  m  or  m  <   0). 
m 


jet   b 

ra,  ,  » . 

^      e  qual 

k 

1_> 

n  n,  -ra.,  n  -m 

.oil  o      q      q 

•    •    •  .A. 


(ra-,  ) ,  » .  (m^)    a^  ^      x, 

1  q        n,  , .  c .  ,  n_      1 


n=0   n,+c..+n  =n        ""  ^        -^i*  •••»-□ 

1  q 

then 

[sk]   ,  ^  m^  ^  m 

1  q 

The  error  in  the  approximation,  i.e.  F(x.)-f  ,  (^°), 
can  be  separated  into  two  kinds.   Part  of  the  error 
comes  from  the  fact  that  the  approximation  is  a  polynomial 
and  that  the  coefficients   of  order  greater  than  [sk] 
are  zero.   This  is  called  the  error  of  the  first  kind 
and  is  denoted  by  E, .   The  rest  of  the  error  comes 
from  the  fact  that  the  coefficients  of  this  polynomial 
approximation  differ  from  the  corresponding  coefficients 
of  the  exact  solution.   This  is  called  the  error  of  the 
second  kind  and  is  denoted  by  Ep .   The  error  obtained 
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when  applying   the   second  method  of   continuation, 
described   in  Section  III,    can  also  be   divided   into 
errors    of  the   first    and   second  kind.      Prom  this    it 
follows   that  Fix^)   =  f^   ^^  ^1  ^'^^l'^^2  ^^®^® 

CD  m.,  in 

o      1  e-o      q 


and 

fek]  ,  m,  m, 

Ep=SZIZ  (b  ^   -b^  ^    )€°      ^...4° 

2        t^^..,^=m     ^i>'''>\     m-L,...,m        1  q 

1  q  ^  ^ 

Since  E,    is    just   the    tail   end  of   the   expansion  for 
f {£,.),    by   assumption,    it  vjill    converge    to    zero. 
Therefore    the   problem  is   to    show  that   £„   ->   0   as   k  ->   od 
Prom  the    assumptions   on  the   domain  of   convergence   of 
P(x. ),    it   follows   that   we   may   use    Cauchy's    Integral 
Formula   over   and   over   to    obtain  the    inequality   on 
the    coefficients: 


la         I  <      " 


q 


n-| ,  .  o « ,  n  '  —   n,      n 

1     q 

where  M  is  the  maximum  value  of  P(x. )  for  |x. I  <  R. , 

see  [10 Jo   Using  this,  it  follows  that  |Ep|  is  less 

than  or  equal  to 

[sk]  CO 

ZIIZ  JZ        ZZ  MCj;i).,.(^q)    A^B- 

m=0  mT+oo.+m  =m  n=k+l  n-,+.,o+n  =n  1  q 

1  q  1  q  ^ 
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It   is  V7ell   known  that 

2=     ^       (hI)--(^'  =  (V-;-%)  .  (1) 

itIt  + .  . .  +m  =ni        1  q 

1  q  ^ 

(This  may  be   found   in  most   texts   on  elementary  probability. ) 
The    following  is   also    true: 

V~  1   =    (^-^q-l)    .  (2) 

n-+. . .+n  =n  ^ 

1  q 

(The  binomial    coefficient  may   easily  be    shown  by 
induction   on  n   to    satisfy  the    same   recursion  relation 
that   the   sum  does.)      Putting    (l)    and    (2)   together  we 
obtain: 

lEpI    <  YZ^       M(^)(^-'2-l)    A^   B™    .  (3) 

'^     ~  m=0   n=k+l       ^  ^ 

Case   one:      It    is  possible   that   g(z)   <  1   for   all 

0   <   z   <   1.      In  this    case    the    approximation  will    converge 

for    any    s,    0   <    s    <   1.      Since    g(z)   <   1,    A(B+1)    <   1. 

But   this   implies   that    Ix?  1  + Ix^'-x?  I    <  R.    for   all   i 

or    Ix.  I    <  R.  ;    that    is,    x.    is   in  D   ,    and  the    case    is 
'    1 '  i'  '       1  o* 

essentially   trivial.      If  we   let    the  m   index  run  from 
0   to   n  we   obtain: 


00  n 

y-        M(^^^-l)    A^'f^    (^)    B^ 


"^  n=k+l  "  m=0 


oo 


<  JZ        M(^''2-^)(A(1+B))^ 
n=k+l  ^ 
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The  expression  on  the  right  is   clearly   the  tail  end 
of  a  power  series  which  converges  if  A(1+B)  <  1, 
Hence,  in  this  case,  Ep  ->  0  as  k  ->  cd  . 

Case  two:   g(z)  =  1  for  some  z,  0  <  z  <  1. 
Since  the  maximum  of  g(z)  is  achieved  at  z  =  B/(1+B), 

s  <  B/(l+B)  .  (1+) 

The  ass-umptlon  about  s  implies   g(s)  <  1.   Prom  the 
inequality  g(s)  <  1  and  (I4.)  It  follows  that 

s  <  1-A  .  (5) 

Rearranging  terms,  we  may  write  (3)  in  the  following 

way : 

00     ^  ^.„   -,       [sk]  ^ 

"^  n=k+l      ^    m=0  "^ 

Consider  the  expression  "5    (  )  B  .   (r  ,  i)  B     will 

^ —r^      m  L  sk  J 

m=0 

be   the    largest   term   in  this    sijtm   since 

/       n       ^   p.[skj-l        /     n     N   ^[sk] 
^[skJ-1^   ^  -   4skJ^   ^ 

°^  [sk]      ^   p, 

n-LskJ+1  -  ^   • 

This   follows   immediately   since,    by    (I4.),    B  >   l/(l-s) 

and  n  >   k+1    implies   l/(l-s)   >    [ sk J/(n-[ skJ+1) , 

Therefore    it    follows   that 

^     "1  -  LskJ 

and 

CO 


IE2I  <  MCsk+DB^^^J^::  (r3ki)(''T^^^''- 


n=k+i    t^^^ 
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Changing  the   index   of   siiramation  this   may  be  rewritten: 

lE^l  <  M(sk+l)B^^(^^l)(k;j)  ^k+1  ^  ^^n 

where 

/ n+k+1 \ , n+k+q  ^ 
w  =   sk   ^^  n+k+1  ^ 
.k+lwk+q^ 
^sk  ^^k+1^ 

[sk]  has  been  replaced  by  skj  this  simplifies  the 
notation  and  It  does  not  change  the  analysis « 

u  =  (k+1    +l).o.(k+l    +n)(k+q+l). »o(k+q+n) 
^        (k+l-sk+lj.,»(k+l-sk+nj(k+l+lj, .. (k+l+n; 

<  (-1-)^  (^13,)^ 

k+n     1      1   k+q+n   k+q 
This  follows,  since   ^.3^^^  <  ^T^     and  j^:^  <  ^ 

for  all  n  >  0.   Hence 

n-0 


But  the  infinite  siom  converges  since  (5)  guarantees 

that   T *  TTT^  <  1  if  k  is  large  enough.   Therefore 

1-s      k+1  °  ^ 


^"    (l-s)(k+l) 

Using   Sterling's   Approximation  to    the   binomial 
coefficient    (      ,  )    ,    we    finally   obtain: 

1^    ,    ^        M(sk+1)       ,k+q^  (g(3))^ 

^2^^^^!-^^  1-    (1-s    (k+1) 
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But  (t^+?)  is  just  a  polynomial  in  k,  therefore  the  whole 
expression  will  be  dominated  by  g(s))  o   Since 
g(s)  <  1,  Ep  ->  0  as  k  ->  00  o 

It  should  be  noted  here  that  E-  converges  to 
zero  under  the  assianptions  that  x     is  in  D  and  that 

^  0  0 

s  satisfies  the  condition  that  s  is  less  than  the 
roots  of  the  equation  g(z)  =  1.   It  is  not  necessary 
for  the  convergence  of  Ep  that  x-,  is  in  D,  o   This 
fact  will  be  used  later  in  showing  that  the  second 
method  of  continuation  convergeso 

Concerning  Repeated  Continuations.   Values  of  an 
analytic  function  obtained  through  repeated  analytic 
continuation  may  be  approximated  by  truncated  power 
serieso  For  simplicity  the  method  will  be  stated  for 
analytic  functions  of  only  one  variable.   Since  the 
proof  follows  closely  that  which  was  already  given, 
just  the  statement  follows o 

Let  f (x)  be  an  analytic  function  given  as  a  series 
in  powers  of  x<,   Let  x-.  ,  o  r.  o  ,x     be  a  sequence  of  points 
about  which  the  power  series  may  be  analytically 
continued,  ioe<,  the  power  series  may  first  be  rewritten 
in  powers  of  x-x, ,  this  new  power  series  may  be 
rewritten  in  powers  of  x-Xp,  etc..   Using  the  sequence 
of  continuations  indicated,  the  problem  is  then  to 
find  an  approximation  to  f(x  )« 
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Call  the  origin  Xq .   Then,  let  R^  be  the  radius  of  conver- 
gence  of  the  power  series  about  Xq,  let  R2  be  the  radius 
of  convergence  of  the  power  series  about  X]_,  etc..   From 
the  assumptions  about  how  the  function  may  be  analytically 
continued,  it  follows  that  |x- -x. _, [  <  Rj ,  ^or  i  =  1  to  q. 
It  is  assximed  that  only  the  terms  of  order  less  than 
or  equal  to  k  of  the  power  series  about  the  origin  are 
known*   Let  s.,  1  =  1  to  q-1,  be  a  sequence  of  real 
nxjmbers,  0  <  s,  <  1.  The  k'th  order  polynomial,  which 
is  an  approximation  to  the  power  series  about  the 
origin,  is  rewritten  in  terms  of  powers  of  x-x, , 
and  terras  of  order  greater  than  (s-k]  are  discarded. 
This  new  polynomial  is  rewritten  in  terms  of  powers 
of  x-Xp,  and  terms  of  order  greater  than  ls-,S2lcJ  are 
discarded,  etc.*   This  process  is  repeated  until  a 
polynomial  of  order  [s,.*.  s   ^kJ  In  powers  of  x-x  , 
is  obtained.   Then  x^  Is  substituted  for  x«   If  the 
s.'s  are  properly  chosen,  the  above  approximation  will 
converge  to  f(x  )  as  k  ->  00 • 

Assuming  the  s.*s  are  given,  define  the  A. 's  by 
Induction: 


Aq  -  1 


A.  =  i . 

1     s.      1-a. 


a^  ^(l-Sj^)   ^ 


*m"^i 


*r*i-i 


S4 

A  'i 
^1+1  • 


-  21  - 


If 


1-s, 


^i+l"^i 


^i-^1-1 


s.  <  1  - 


1  1-1 


and  if  A, 


1   1-1 


R. 

1 


<  1   for  1  =  l,oo«,q-l, 
then  the  approximation  converges  as  k  ->  oo  « 
Purthermore,  given  any  set  of  x. 's,  and  R. 's,  there 
exist  s.'s  which  satisfy  the  above  conditions. 
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3.  The  second  Method  of  Continuation, 

Description  of  the  Method,   This  method  makes  more  use 
of  the  differential  equations  than  the  first.   In  order 
to  simplify  notation  the  method  will  be  described  and 
the  proof  of  convergence  will  be  given  for  an  analytic 
initial  value  problem  in  only  two  variables  and  one 
unknown.   With  only  trivial  modifications  the  method 
and  proof  are  applicable  to  the  most  general  analytic 
non-characteristic  initial  value  problem. 

and  u(0,y)  =  f(y),  f(0)  =  0  (6b) 

OP 

where  G(u)  is  analytic  at  u  =  0,  G(u)  =  >    a  u  , 

n=0  00^  -  ^ 
and  f (y)  is  analytic  at  y  =  0,      f (y)  =   >   by   , 

ri=0 

and  only  the  coefficients  of  order  less  than  or  equal 
to  k  of  this  expansion  of  f (y)  are  known.   Let  u,  (x,y) 
be  the  polynomial  which  Is  obtained  from  u(x,y)  when 
u(x,y)  is  expanded  in  a  series  in  powers  of  x  and  y 
and  terms  of  order  greater  than  k  are  discarded. 

Prom  the  differential  equation  (6a)  it  is  possible 
to  write  the  partial  derivatives  of  u  of  any  order  in 
terms  of  partial  derivatives  of  u  with  respect  to  y 
of  that  order  or  less  by  differentiating  over  and 
over  again.   Therefore,  using  (6b),  it  is  possible  to 
compute  u,  (x,y)  when  only  the  coefficients  of  f  (y) 
of  order  less  than  or  equal  to  k  are  known. 
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Let    (x    ,Y    )  be    a  Doint    in  the   domain  of   absolute 
o      o 

convergence    of  the   power   series    of  u(x, y)    expanded 

about    the   origin   (the    choice    of   this    point  will  be 

further  restricted  to   insure    convergence   of  the  method  ) 

and  let    (x-,,y-,)   be    a  point    in  the    domain  of  absolute 

convergence    of   the   power   series    of  u(x,y)    expanded   about 

the   point    (x    ,y    ).      The   problem    then   is    to   find   an 
o      o 

approximation  to   u(x-,,y-,). 

Introduce   the    new    coordinate    system   E,   =  x-x 

and    ^  =  y-v    ,    and   let    ■£      =  x, -x      and   ^     =   y^ -y    • 
^        -^"o'  o  lo  ^o        "^l-'o 

Let  ^  =   G(w)  ^  (7a) 

and  w{0,K)    =  u^(x^,^+y^)  (7b) 

and  finally,  let  Wr  ,1(4,^)  be  the  polynomial  which 

is  obtained  from  w(^,^)  when  w(4,^)  is  expanded  in 

powers  of  ^   and  ^   and  all  powers  of  order  greater  than 

[sk]  are  discarded  where  s  is  some  constant  between 

zero  and  one.   ([sk]  is  the  greatest  integer  in  sk. ) 

If  equation  (7b)  were  replaced  by  w(0,^)  =  u(x  , ^+y  ) 

then  w(^  ,^  )  would  equal  u(x, ,y, ).   As  it  stands, 

we  shall  show  that  Wr  ,  i(£  ,^  )  ->  u(x-,  .y^  )  as  k  ->  00 

[skj  ^0*^0        1**^1 

if  s  and  (x  ,y  )  satisfy  certain  restrictions. 

These  restrictions  on  s  and  (x  .y  )  can  be  stated 

o'  o 

in  terms  of  the  following  dominating  problem: 
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%  -   m.)  ^  (8a) 


and  v(0,7)  =  f(7)  (8b) 

_  00  --  _  00 

where   G(v)  =  J~~    \&    \    v   and  fij)    =  2ZI  I^  \j      • 
n=0   ^  n=0   ^ 

The  terra  "dominating  problem"  requires  some  explanation. 

Let  ^{^,K)=   Z  Cin,n^''^  ^^  ^^^'^^  =  ^  ^m,n^'^^''-  ^^^'^^ 
is  said  to  dominate  ^i{^,K),    (I)(4,^)^  t(?,^),  if 

c    >  Id    1.   A  problem  (a  differential  equation  plus 
m,n  —  '  m,n'     ^  ^       i- 

an  initial  condition  in  the  same  form  as  (6a)  and  (6b)) 

is  said  to  dominate  another  problem  if  the  functions 

occuring  in  the  former  dominate  the  corresponding  functions 

occuring  in  the  latter.   This  concept  is  useful  for  it 

is  easy  to  show  that  the  solution  of  the  dominating 

problem  dominates  the  solution  of  the  other  problem. 

Since  G(u)  ^  G(u)  and  f (y)  D  f(y),  problem  (8)  dominates 

problem  (6),  hence  v(x,y)  ^  u(x,y).   This  concept  of 

domination  is  the  most  important  tool  in  the  proof  of 

convergence. 

Theorem:   If  (x  ,y  )  is  in  the  domain  of  convergence 

of  the  power  series  of  v(x,y)  expanded  about  the  origin 

and  if  (x,  ,y-.  )  is  in  the  domain  of  convergence  of  the  ^ 

power   series    of  u(x,y)    expanded   about   the    point    (x    ,y    ) 

then   there   exists    an  s    such  that  Wr    i^  i  ( C   tK    )    ~>  '"-i(x-|  >7-]  ) 

as   k   ->   00  . 

Since  (x  ,y  )  is  in  the  domain  of  convergence 
o   o 

-25  - 


of  the  power  series  of  v(x,y)  expanded  about  the  origin, 
there  exist  two  positive  numbers,  X  and  Y,  such  that 
(X,Y)  Is  In  this  domain  and  |x^!  <  X  and  ly  1  <  Y. 

Let  A  =  max(— ^T — ,—y — j,  and  B  =  max{-r^ — ['"\'y   'I  ' 

and  g(z)  =  AB  — =j — -  , 

°  z,-,       \l-z 

z  (1-z) 

If  s  is  chosen  such  that  It  is  positive  and  less  than 

the  roots  of  g(z)  =  1,  then  the  approximation  will 

converge o 

The  domain  of  absolute  convergence  of  the  power 

series  of  v(x,y)  expanded  about  the  origin  may  well  be 

properly  included  in  the  corresponding  domain  of  u(x,y)o 

The  restricting  of  (x  ,y  )  to  this  possibly  smaller 

domain  is  unnecessary  in  the  proof  of  the  convergence 

of  the  first  method  of  continuation  and  seems  somewhat 

arbitrary  in  the  present  caseo   (An  example  has  been 

found  for  which  this  restriction  is  unnecessaryo ) 

However,  since  the  Cauchy-Kowalewski  Theorem  guarantees 

that  this  possibly  smaller  domain  is  not  empty,  the 

second  method  of  continuation  will  converge  for  a  larger 

domain  than  the  domain  of  absolute  convergence  of  the 

power  series  of  u(xs,y)  expanded  about  the  origin  unless 

the  boundary  of  this  domain  is  the  natural  boundary  of 

u(x,y)o   That  is,  even  though  the  second  method  of 
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continuation  can  only  be  shown  to  converge  in  a  domain 
which  may  be  smaller  than  the  domain  of  convergence  of 
the  first  method,  this  domain  always  contains  points 
for  which  the  functions  could  not  be  approximated  unless 
some  method  of  continuation  were  used. 
Proof  of  Convergence.   It  is  convenient  to  introduce 
five  new  operators,  S,  I,  I,  T,  and  R,  and  then  to 
describe  the  approximation  in  terms  of  these  operators. 
The  operation  of  substituting  x  for  the  variable 

X 

X  in  a  function  (|)(x)  is  indicated  by  S  °(|)(x);  that  is, 
s/(l)(x)  =  (!)(x^). 

The  operator  I  (l  for  initial  value  problem)  may 
only  be  applied  to  a  function  of  ^.   By  I(j)(^)  is  meant 
the  unique  analytic  function  of  the  variables  ^   and  ^ 
which  solves  the  partial  differential  equation  (7a) 
and  is  equal  to  (f)(^)  when  ?  =  0.   That  is,  \|/(?,^)  =  I(t)(^) 
if  and  only  if  ^  =  G(\|/)  tt  and  \l/(0,^)  =  <^{K)'      The 
operator  I  is  the  same  as  I  except  the  differential 
equation  -jp  =  G(\l/)  4^  is  used  Instead  of  (7a). 

The  operator  T.*'^  (T  for  truncate)  can  only  be 
applied  to  an  analytic  function,  ^ix,j),    which  is  analytic 
at  the  point  x  =  0,  y  =  0.   By  T'^.^'^^ix^j)    is  meant  the 
polynomial  in  x  and  y  which  is  obtained  when  (|)(x,y)  is 
expanded  in  powers  of  x  and  y  and  terms  of  order  greater 
than  j  are  discarded.   That  is,  if 
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q=0  rri+n=q        * 


in     n 


then 

T^*^(l)(x,7)   =2^X1        ^^n^"^^ 
J  q=0  m+n=q     ^'^ 

T.*^.  T.         ,  etc.  are  defined  in  a  similar  manner, 
J   *   J 

The  operator  R  (R  for  remainder)  may  be  applied 
to  the  same  functions  that  T  can.   R''!"*^(f)(x,y )  is 

J 

defined  to  be  equal  to  <^{x,j)    -   T^»^(|)(x,y ) ,  etc.. 

Let  u(x,y)  be  the  solution  of  (6)  and  let 

U(4,^)  =  u(C+x  ,^+y  ),   Let  v(x,y)  be  the  solution 

of  (8)  and  let  V(C,^)  =  v( ?+ |x^ 1 , ^+Iy^ | ) . 

£+!x  1   ^+1y  I 
Let   vi(4,^)  =S^   °'(S^   °  (Tj^^v(x,y))), 

UgC?,^)  =  I(S^°(S^'^°(T^^yu(x,y))))  . 

Sg.°S  °T^*^,U2(€,^)  is  the  approximation  already 
described.   The  problem  is  then  to  show  that 

as  k  ->  oo .   u(x-,,y,  )  can  be  written  in  the  follovjing  way; 
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^      K      r    y 

Since  (x-,j,y^)  was  assuiried  to  be  in  the  domain  of 

convergence  of  the  power  series  of  u(x,y)  expanded 

about  the  point  (x  j,y  ),  this  implies  that  (^  ,  ^  )  is 

in  the  domain  of  convergence  of  the  power  series  of 

4  ^  €  ^ 
U(4,^)  expanded  in  powers  of  C  and  ^o   S^  S^  ^[skl^^^»^^» 

the  error  of  the  first  kind  as  described  in  Section  II, 


is  just  the  tall  end  of  the  power  series  for  U(^  ^^  ), 
and  therefore  by  assumption  converges  to  zero.   Therefore 
the  problem  is  reduced  to  showing  that  the  error  of  the 
second  kind  converges  to  zero,  i.e. 

S  °S  °T|-*^j(U(C,^)  -  yx^{^,K)) ->   0  as  k  ->  oo. 

This  is  easily  established  by  the  use  of  two  simple 
lemmas. 

Lemma  1:   V(4,^)  ~  Vgl^,^)^  U(4,^)  -  u^C?,^). 

U  !    ^+iy  I 

Let  M(^,5)  =  S^  °  S^    °  (Tj*^v(x,y)  +  ml»^vU,j)) 

X  ^+y 
and  N(^,5)  =  S^°Sy   °(Tj*^u(x,y)  +  5R^»^u(x,y ) ) . 

Since  v(x,y)  ^  u(x,y),  it  follows  that 

M(^,5)  3  N(^,5)  .  (9) 
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The  differential  equations  associated  with  I  and  I 
can  be  thought  of  as  containing  5  vacuously;  therefore, 
T(M(^,5))  ^  I(N(^,5))  with  respect  to  the  variables 
C,  K,    and  5.   Since  3°  l(M(c;,5))  =  v^{^,t,)    and 
Sg  I(n(^, 5))  =  UpCC,^),  there  exists  two  analytic 
functions  M^(?,^,5)  and  N^(?,^,5)  such  that 

T  (M(^,5)  =  v^UyK)    +   5M^(C,^,5)  (10) 

and 

I(N(^,5))  =  u^i^,K)    +  5N^(C,2;,5)  (11) 

Prom  (9)  it  follows  that  M^(^,^,5)  2   N^(^,^,6).   This 
would  still  be  true  of  1  (one)  were  substituted  for  5. 
Substituting  1  for  5  in  equations  (10)  and  (11)  we 
obtain: 

V(C,^)  =  V2(4,^)  +  14-^(^,^,1) 

U(^,^)  =  u^i^,V    +  N^(4,^,l) 

Therefore      V(C,^)    -  v^{^,K)  2U(?,^)    -  n^ii,K)» 
Lemma   2:      V(C,^)  2  V2(?,^)   2  V3_(?,^)  2  ^    • 

Since   v(lxj,^+!yj)  ^  \  °  S^        °  Tj'yv(x,y), 
it   follows   that 

T(v(lxJ,    ^H-IyJ))2T(S^''°'s^"''^°'Tj'7v(x,y)). 

But        V(?,^)   =  T(v(Ix^l,^+Iy^I)) 

and        V2(€,^)   =  Ks^  °   S^        °  T^'^v(x,y)), 
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hence  V(?,^)  ^  v^C?,^)  . 

=  t^'^''^°'''^o't(t,^'''^°'''^°'w|xJ,,.|.J))       (12) 
=  T^'''''°''''°'l(v(|xJ,5.|yJ))  (13) 

I 


=  S^   °   Sy   °  T^^y  v(x,y) 

Vp(C,^)  is  an  analytic  function  In  \x    \    and  ly  |  as 
well  as  E,   and  ^,  therefore  it  is  meaningful  to  truncate 
v^CC,^)  with  respect  to  <g,  ^,  Ix^|,  and  ly^U   In  going 
from  (12)  to  (13)  the  truncation  of  v( |x^ 1 , ^+ jy^ I )  is 
dropped.   Since  only  the  coefficients  of  v  of  order 
greater  than  k  are  disturbed  in  the  initial  conditions 
of  the  I   problem  in. (12),  only  the  coefficients  greater 
than  k  of  the  solution  of  the  I  problem  are  disturbed. 
Since  the  solution  Is  truncated,  the  original  truncation 
has  no   effect.   Since 

^,'is  U^l,  ly,  ■ 


T,_     ^°  '      °   V2(4,^)  =  v^i^^K) 


■k 


and  Vp(^,^)  has  only  positive  coefficients  considered 
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as  a  function  of  ?,  ^,  I^o'*  ^^^  '^o'*  ^^   follows  that 
Vp^S*^)  Z  "^l^^*^-^"   v^(C,^)  has  only  positive  coefficients, 
therefore  it  dominates  zeroo 

Putting  Lemmas  1  and  2  together: 

Therefore 

Is  greater  than  or  equal  to 

s^^°'s^^°'t^;^^  (v(?,^)  -  v^i^,V)  (15) 


is  greater  than  or  equal  to 

(II].)  is  the  error  of  the  second  kind  when  the  first 
method  of  continuation  is  applied  to  the  dominating 
problemo   (15)  is  the  error  of  the  second  kind  when 
the  second  method  of  continuation  is  applied  to  the 
dominating  problem,   (16)  is  the  absolute  value  of  the 
error  of  the  second  kind  when  the  second  method  of 
continuation  is  applied  to  the  original  problem.   Since 
(1^  U ly  1)  i2  i^  the  domain  of  convergence  of  v(x,y) 
expanded  about  the  origin  and  s  is  appropriately  chosen, 
(iLj.)  converges  to  zero  by  Section  IIo   Therefore  (l6) 
converges  to  zero,  which  completes  the  proof,, 
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(16) 


With  almost  no  changes  in  the  statement  and  proof 
of  convergence,  the  second  method  can  be  expanded  to 
include  any  system  of  quasi-linear  partial  differential 
equations  in  any  number  of  -unlcnowns.   It  should  be 
pointed  out  that  this  also  includes  ordinary  differential 
equations.   The  domain  of  convergence  of  the  power  series 
of  one  unknown  about  some  point  must  be  altered  to  mean 
the  largest  open  set  in  which  all  of  the  power  series  of 
the  unknowns  expanded  about  this  point  converge  absolutely, 
Also,  the  definition  of  A  and  B  must  be  appropriately 
altered. 

Since  (ll|)  is  greater  than  or  equal  to  (15)  one 
might  suspect  that  the  second  method  of  continuation 
is  better  than  the  first »   This  actually  proves  to  be 
true  in  some  simple  examples  as  well   as  in  the  computa- 
tions done  on  the  detached  shock  problem.   These  results 
ar'e  described  in  the  next  sectlon-o 
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I|-»  Results  of  the  Test  Calculations. 

Description  of  the  Detached  Shock  Problem.   Let   z,r 
be  cylindrical  coordinates.   It  was  assuitied  that  the 
flow  is  symmetric  about  the  z-axis,  and  that  the  shock  is 
given  by  z  =  H(r)  where  H(r)  is  an  analytic  function. 
In  the  region  ahead  of  the  shock,  given  by  z  <  H(r),  the 
flow  quantities  are  constant  and  the  flow  is  parallel 
to  the  z-axis.   The  shock  was  ass-umed  to  be  a  hyperboloid 
of  revolution,  the  gas  was  assumed  to  be  polytropic 
with  a  Y  of  '^"k-f    ^-^^^  M  (the  mach  number)  was  taken  to 
be  12,   Since  the  publishing  of  the  earlier  report 
(6),  it  has  been  pointed  out  by  G«  B,  Whitham  and 
others  that,  with  such  a  large  Mach  number,  y  ~  1«U- 
is  not  realistic,  owing  to  dissociation  of  the  gas. 
In  order  to  be  realistic  y  ~  ^'2.   would  be  a  better 
approximation  for  air.   However,  since  the  primary 
interest  in  the  computation  was  theoretical,  y  ~  ^'h 
was  retained  so  that  the  present  results  could  be 
compared  with  those  from  the  Univac. 

In  order  to  apply  the  methods  described  it  is 
necessary  to  introduce  a  curvilinear  coordinate  system 
which  has  the  property  that  one  coordinate  is  constant 
along  the  shock.   This  is  the  only  property  that  the 
coordinates  must  satisfy  other  than  that  the  functions 
involved  in  defining  the  coordinate  system  must  be 
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analytic.   In  order  to  achieve  this  the  coordinates 
X  and  y  were  introduced: 

X  =  z-H{r)  z  =  x+h(y) 

or 

■^  o  "^      o 

where    r      Is    a    constant   and  h(y)    =   H(y+r    )„      These 

particular    coordinates    have    the    property   y   is    a  radial 

coordinate   measured   from   the   point    r   =  r       (not   from   the 

o 

axis)  and  x  is  a  coordinate  measured  parallel  to  the 
axis  from  the  shock  (not  from  a  fixed  plane).   The  flow 
quantities  were  so  scaled  that  the  radius  of  ciorvature 
at  the  tip  of  the  shock  was  1.   This  guarantees  that  the 
tip  of  the  body  is  located  at  z  =  ,1,  r  =  0  approximately. 

Let  p  be  the  density,  p  the  pressure,  u  the 
component  of  velocity  parallel  to  the  axis  of  symmetry, 
and  V  the  component  perpendicular  to  the  axis.   Using 
the  Rankine-Hugonlot  shock  conditions  the  flow  quantities 
(u,  V,  p,  and  p)  can  be  determined  as  functions  of  y 
when  X  =  0,   The  differential  equations  are  easily  derived 
from  the  law  of  conservation  of  mass,  the  law  of 
conservation  of  momentum,  the  assiomption  that  changes 
are  adlabatic  (that  is,  entropy  is  constant  along  stream- 
lines), and  the  equation  of  state  for  a  polytroplc  gas. 

Most  of  the  computations  were  directed  toward  finding 
the  location  of  the  body.   The  surface  of  the  body  can  be 
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defined   in  terms    of  the    stream  function,    \|/«      t  may  be 

determined   in  terms   of   the   other   flow  quantities  by 

means   of   quadratureso      In  this   problem  the    streamline 

\1/  =   0   presumably   consists    of  two  parts:      the    z-axis 

and   a   cixrve    intersecting  the    z-axis    at   a  positive   value 

of   Zo      The    second  part  may  be   taken  as    defining  the 

surface    of    the   bodyo 

r     was    chosen  and   the   flow   quantities    (u,    v,    p, 

and   p)   were    approximated  by  truncated  power   series   in 

X   and  y*       (in  terms    of  the    z,r   coordinate    system  the 

expansions  were  made    about   the   point    (H(r    ),r    ); 
•^  ^  o        o 

changing  r      allows   one    to    investigate   different  parts 
of   the    surface    of   the   body^ )      Using   the   notation  of 
Section  III,    (x   ^y    )   was    chosen:      and,    then,    using 

CO 

either  the  first  or  the  second  method  of  continuation, 
the  flow  quantities  were  approximated  by  truncated 
power  series  in  ^   and  ^  where  £,   =   x-x   and  ^  =  J'J^" 
Both  the  truncated  series  in  x  and  y  and  in  E,   and  ^  were 
computed  to  a  fixed  order  k«   In  most  of  the  computations 
k  was  taken  to  be  19  (the  largest  possible  k  due  to 
machine  restrictions ) «   The  truncated  power  series  in 
<C  and  ^  of  the  flow  quantities  were  then  further  truncated 
so  that  all  of  the  terms  of  order  greater  than  [skj  were 
discarded.   Using  these  approximations  to  u,  v,  p,  and  p, 
an  approximation  to  ^(?,^)  was  computed^  call  this 
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^^{E,,^).      ^     was    chosen   and    finally   ^    ,    the    solution  to 
\|;-(^,^    )    =   0,    was    cornputedo      In  the    original   z,r 
coordinate    system,    then, 

(V^o-^-^^^o-^^o^^yo-^V^o^ 

is  an  approximation  to  a  point  on  the  surface  of  the 

body» 

Results  Using  the  Second  Method.   In  order  to  compute  a 

point  on  the  body  surface,  six  constants  (k,  r  ,  x  ,  y  , 

[sk],  and  ^  )  must  be  chosen.   This  number  was  too  great 

to  allow  a  thorough  investigation  of  the  effect  on  the 

accuracy  of  varying  these  constants.   In  all  of  the 

results  discussed,  k  =  19.   The  constants  for  which  the 

best  results  were  obtained  are  indicated  schematically 

in  Figure  3  page  39.   For  example,  with  the  choice  of 

constants  r   =  .Ll,  x  =  o03,  y  =  0,  ^  =  -ol,  first 
o     'o      '•'o       o 

the  truncated  power  series  in  x  and  y  about  the  point  A 

were  obtained,  then  truncated  power  series  in  £,   and  t, 

about  the  point  B,  and  then  the  body  point  C,   ^  was 

always  chosen  a  multiple  of  o05,  and  y  was  always 

taken  to  be  zero.   The  best  results  were  obtained  by 

the  following  combinations  of  x   and  r  :  r  -    o2, 
^  o      o    o 

^o  "  '^^'^    ^o   "  °^'    "^o   "  '^^'^   ^o  "  '^*  ^o  "  '°^'' 

r  =  e8,  X  =  ,08.   Although  few  cases  were  tried  other 
o     '   o  '^ 

than  these  alreadv  listed,  it  seemed  clear  that  (x  ,y  ) 

3  o   o 
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Fig-ure    3 
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must  be  taken  closer  to  the  shock  than  to  the  bodyo   The 
point  E  in  Plg-ure  3  requires  some  explanation.   The 
truncated  power  series  about  this  point  was  not  obtained 
directly  from  the  point  on  the  shock,  but  was  obtained 
through  the  use  of  the  second  method  of  continuation 
from  the  truncated  power  series  about  Do   The  results 
were  never  improved  by  taking  [skj  less  than  19 „ 

In  choosing  the  best  result,  three  checks  were 
available : 

1)  Bernoulli's  law  states  that  (u^+v^)/2  +  yp/((y-1)p) 
should  be  constant  (208o6  for  this  particular  problem) 
throughout  the  flow„ 

2)  Consider  some  point  in  the  flow  where  the  pressure 
and  density  are  p  and  po   Let  p..  and  p,  be  the  pressure 
and  density  of  the  flow  just  after  the  streamline  which 
goes  through  this  point  crosses  the  shock.   Since  entropy 
is  constant  along  streamlines,  p/p-,  should  be  equal  to 
(p/p^ )  .   The  details  of  computing  p,  and  p^  are  described 
in  Appendix  Ao   This  check  proves  to  be  a  little  too 
sensitive,  probably  due  to  the  great  number  of  approxima- 
tions involved. 

3)  The  flow  quantities  at  a  point  may  be  computed  with 

different  choices  of  r  ,  x  ,  etco ;  these  of  course  should 

o    o 

agree o 
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To  illustrate  these  checks  the  following  example 
is  given: 

r       z        Ber,       (p/p-j^)^      p/p-^ 

(1)  .05   .10328060  20806OOOO  1.1015690   1.1015690 

(2)  .05   .10328060  208.60001  1.1015696   1.1015695 

The  entrances  in  row  (l)  were  computed  with  r   =  .2, 

X  =  «03,  and  Z,     =   -.15;  the  entrances  in  row  (2)  with 
00       ' 

r  =  0,  X  =  .03,  and  ^  =  .05.  7     =  0  in  both  cases. 
o     'o  o  '^O 

It  is  not  unreasonable  to  assume  that  the  figin'es  given 
in  column  z  have  eight  digits  of  accuracy.   Of  coTirse, 
there  is  always  the  possibility  that  these  figures  are 
deceiving.   But  this  could  always  be  true  for  any  method 
for  v;hich  there  is  no  exact  check  on  the  answers.   In 
the  present  case  it  seems  very  unlikely. 

Table  1  on  page  62  gives  the  location  of  the 
surface  of  the  body  in  the  z,r  coordinate  system  for 
various  values  of  r  and  also  the  flow  quantities.   Only 
the  digits  which  were  felt  to  be  accurate  are  listed. 

It  should  be  pointed  out  that  the  power  series 
which  were  computed  bear  out  the  conjectiore  that  some 
method  of  continuation  must  be  used  in  order  to  compute 
the  surface  of  the  body.   When  r  =0,  the  first  20 
coefficients  of  u(x,0),  p(x,0),  and  p(x,0)  behave 
approximately  like  the  first  20  coefficients   of 
l/(.08+x).   Therefore  the  body  would  seem  to  be  outside 
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the  domain  of  convergence  of  the  power  series  of 

u(x^y)e   This  same  property  seems  to  be  true  for  all 

the  truncated,  power  series  expanded  about  points  on  the 

shock  surface. 

Results  of  the  First  Method,   Since  the  results  using 

the  second  method  seem  to  be  accurate  to  eight  digits, 

these  results  were  used  as  a  guage  for  checking  the 

accuracy  of  the  results  using  the  first  methodo   All 

of  the  computations  using  the  first  method  used  k  =  19 

and  r   =  o2.   Various  choices  for  x   and  y  were  used: 
o  o     •'o 


«03,  y^  =  -.1;  x^  =  .03,  y^  =  o;   x^ 


■05,  y^  =  0; 


X     =o03,    y     =olo      In  all   cases  which  were    computed, 

0  0 


[skj  =  ; 
results, 


,  x^  =  o03,  and  y^ 


,1  gave  the  most  accurate 


r 

ol 

.2 
.3 


1st  method 

z 
,  10852  2)|0 
,12960302 
,l6509k92 


2nd  method 

z 
.10852381^ 
,12960099 
c  16509961; 


The  first  method,  then,  gives  approximately  6  digits 
of  accuracy. 
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The  Choice  of  x   and  s  in  the  First  Methodo   The  first 
o  — 

method  does  not  involve  in  any  essential  way  the  use 
of  the  differential  equations.   It  can  be  thought  of 
as  a  method  for  "analytically"  continuing  a  power  series. 
If  the  method  were  applied  to  an  analytic  function  of 
one  variable,  in  most  cases  x,  and  k  would  be  giveno 
The  questions  is,  then,  how  to  choose  x   and  s  so  as 
to  minimize  the  error  of  the  approximation.   For 
simplicity  assume  that  the  analytic  function  is  given 
as  a  power  series  about  the  origin,  that  the  radius  of 
convergence  of  this  power  series  is  1,  and  that  the  only 
singularity  that  the  function  possesses  is  at  x  =  ao 
(From  this  it  follows  that  jaj  =  1« )   In  order  to  gain 
values  for  x   and  s  two  important  assumptions  are  madeo 
First,  it  is  assumed  that  there  is  no  systematic 
cancellation  between  the  errors  of  the  first  and  second 
kind  as  described  in  Section  II „   In  that  section  an 
upper  bound  was  derived  for  |Epl  (the  error  of  the 
second  kind).   Although  it  is  not  clear  at  first  glance, 
all  of  the  estimates  used  in  deriving  this  bound  are 
extremely  tight.   The  second  assumption  is  that  the  error 
of  the  second  kind  actually  behaves  like  this  upper  bo-und< 
These  two  assumptions  may  well  always  holdo   Under  these 
assumptions  it  must  be  true  that  the  best  choice  for  x 
and  s  must  be  such  as  to  make  the  exponential  part  of  the 
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error  of  the  second  kind  equal  to  the  exponential  part 
of  the  error  of  the  first  kind: 


X. 


X  -X-, 
o   1 


X 


1 


!  (1-SJ 


X,  -X 

1   o 


a-x 


o 


(17) 


The  best    choice   for  x      and   s  would  then  be   the   value    of 

o 

X  and  of  s  that  minimizes  the  right  hand  side  (or  the 
o 

left  hand  side)  of  (17).   For  any  e  >  0,  the  error  in 
the  approximation  for  large  k  would  be  less  than 

^  ^o 

where  M  is  some  constant  which  depends  on  e  and  f(x). 

Even  if  the  two  assumptions  were  dropped,  the  previous 

statement  would  still  be  true,  but  It  might  be  possible 

to  choose  X   and  s  so  as  to  obtain  an  even  smaller  bound. 
o 

The  values,  x   and  s,  which  minimize  one  side  of  (17)  ape 
given  in  Figure  I4.  on  the  next  page  as  functions  of  x. 


for  various  values  of  ao 


E  = 


Xt  -X 

1  o 


a-x 


is  also  listed  as  a  function  of  x-,  .   In  Figure  l\.   the 
graphs  numbered  1  correspond  to  a  =  -1,  those  numbered 
2  correspond  to  a  =  +i,  and  those  numbered  3  correspond 
to  a  =  1/2  +73/2  i.   In  Tables  2  and  3  on  pages  63 
and  65,  these  values  of  x   and  s  are  used  to  compute  the 
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values  of  various  elementary  functions  at  x,  =  08  and 

X,  =  1.2o   The  machine  also  computed  the  best  choice 

for  [sk]  when  fix),  x  j,  x^  ,  and  k  were  given.   This 

best  choice  differed  so  little  from  the  [sk]  gained 

from  the  fixed  s  computed  through  (I7)  that  it  is  not 

listed.   This  supports  the  belief  that  the  two  ass-umptions 

used  to  derive  (l?)  were  very  reasonable^   The  error  was 

at  times  a  factor  of  10  better  when  the  best  choice  for 

[skj  was  used,  but  this  was  never  consistent.   As  can  be 

seen  from  the  tables,  the  error  seems  to  be  affected  as 

much  by  large  coefficients  in  the  original  pox^rer  series 

(eogo  (l/(l+x)-^)  as  by  a  near  singularity  (e.g.  l/(l+x  )). 

Choice  of  x   and  s  in  the  Second  Method,   Hov;  to  make  the 
o  — —   

best  choice  for  x   and  s  in  the  second  method  is  a  much 

o 

more  complicated  question„   It  seems  to  depend  not  only 

on  the  solution  but  also  on  the  differential  equations 

involved.   In  the  detached  shock  computations  it  seemed 

clear  that  the  best  choice  for  s  was  s  =  1,   It  is  easy 

to  find  examples  in  which  s  =  1^  if  not  the  best  choice 

for  s,  is  the  only  reasonable  choice  for  s  given  a 

particular   k„    x    .    etco       (The   best    choice   for    s    can  be 

o 

unreasonable  in  the  following  sense:   the  use  of  the 
best  [skJ  in  the  truncation  produces  the  smallest  error, 
but  the  use  of  [skJ+1  produces  a  very  large  error.) 
One  of  these  examples  is  as  follows: 
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~  =  =u^    where   u(0)  =1. 
dx 

The  solution  is  u  =  l/(l+x)<,   If  the  second  method  of 
continuation  is  used  for  x-,  =  1,  x  =  /2  -  1  and  s  =  1 
are  the  most  reasonable  choices »   The  error  is  then 
asymptotic  to  (-1)^(1- ^)(  /2  -D^'^-'-o   If  the  first 
method  is  used  then  x  =  1/3  and  s  =  l/3  are  the  best 
choices  according  to  the  last  section  and  the  error  is 
bounded  by  M(-^  +  e  )^  or  M(o791|+8)^ 
for  e  >  0  and  M  large  enough.   This  bound  was  confirmed 
by  actual  test  computationso   Hence,  the  second  method 
is  much  superior  to  the  first  in  this  example » 

Since  in  both  the  previous  example  and  in  the 
detached  shock  computations  the  best  choice  for  s  using 
the  second  method  of  continuation  is  apparently  s  =  1, 
it  should  be  pointed  out  that  there  exist  cases  where 
s  cannot  be  so  choseno   If  the  partial  differential 
equations  are  linear  with  constant  coefficients  then 
the  first  and  second  methods  of  continuation  are  identical, 
Therefore  if  there  is  a  singularity  in  the  initial 
conditions,  there  exist  points  for  which  neither  method 
will  converge  unless  s  is  taken  less  than  oneo 


he 


Concluding   Remarks o      Even  though  the  methods   described   give 
great    accuracy^    in  actual  practice   this    Is  balanced  by  the 
fact   that    a   large    amount    of    computing    time    Is   required,, 
In  order    to    compute    the    truncated  power    series    of    the 
flow   quantities    from   the    differential   equations    and  the 
initial   conditions   when  k   -   19,,    20  minutes   is  required 
on    the    IBM   70L|..      If   all    that    is   required    is    the    stand-off 
distance    of   the    shock  from   the   body^    using    either  method, 
the    time   required   for   the   rest    of   the    computation   is 
negligible o      Therefore    the   first  method  requires    approxi- 
mately  20   minutes    and   the    second  method   I4.O   minutes.      If 
the   arithmetic    (see   Appendix    B)  were    simplified   and   the 
code   rewritten^,    the    time    could  be   reduced  by   at  most    a 
factor   of  2» 

Other  methods    are    known  for    solving  the    detached 
shock  problemo      P.    R«    Garabedian   and  Ho    M,    Lelberstein, 
and  M,    D,    Van  Dyke,    for   example,    have   described   successful 
finite-difference   methods,    see    [11]    and    [8J,       Van  Dyke's 
method   gives   the    stand-=off   distance   to    two   figures   of 
accuracy    (which   is   all   that   Is   required  for   practical 
purposes)    in  5  minutes    on  a  much   slower  machine »      Even 
if  the   two   methods      of   continuation  were    streamlined 
as  much  as   possible,    it    seems    doubtful    if   they   could 
ever   compete  with  Van  Dyke's   computing  time. 

It    is   not   remarkable    that   the    first   method    is 
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insensitive  to  whether  the  equations  involved  are 

elliptic  or  hyperbolic,  but  the  second  method  seems 

Insensitive  also.   In  other  methods,  such  as  Van  Dyke's 

special  consideration  must  be  given  to  the  sonic  line. 

With  the  present  methods  no  effect  whatever  was  noticed 

when  crossing  the  sonic  line. 

Even  though  convergence  has  been  shown,  there  are 

many  questions  which  need  to  be  answered:   Under  what 

conditions  is  the  second  method  preferable  to  the  first? 

Are  the  restrictions  on  s  in  Section  II  the  weakest 

possible  under  the  given  assumptions?   (This  seems  very 

likely.)   Is  the  use  of  the  dominating  problem  in 

Section  III  ever  necess^jry,  and  for  what  conditions  may 

s  =  1?   Kay  the  second  method  be  applied  repeatedly?   How 

should  X  -y   and  s  be  chosen  in  the  second  method  so  as 
o '  "^  o 

to  assure  the  maximxmi  rate  of  convergence? 

Changing  one's  point  of  view  might  also  yield 
interesting  results.   For  example,  is  it  possible  to 
obtain  a  useful  statement  of  the  first  method  in  terms 
of  information  theory?   The  first  method  could  be 
thought  of  as  a  means  for  recovering  an  encoded  message 
(the  power  series  of  a  function)  when  there  is  noise 
(truncation).   This  is  made  plausible  by  the  fact  that 
log  (g(z))  contains  a  term  suggestive  of  entropy  or 
information  content. 
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Appendix  Ao   Details  of  the  Test  Problein» 

The  Equations o   Let  z,v   be  cylindrical  coordinates, 

and  let  the  shock  be  given  by  z  =  H(r)o   It  is  convenient 

to  take  the  units  of  length,  density,  and  pressure 

(these  are  dimensionly  independent)  to  be  the  radius  of 

curvature  of  the  shock  at  the  axis,  the  density  ahead, 

and  the  pressure  ahead.   Then  the  conditions  ahead  may 

be  written  as  p  =  1 

p  =  1 

u  =  U 

V  =  0 
Here  u  and  v  are  axial  and  radial  components  of  fluid 
velocity  and  U  is  a  constant  related  to  the  Mach  number 
by  the  equation  v 

U  =  mVy 
where  y  is  the  constant  of  the  polytropic  gaso 

It  was  assumed  that  H(r)  satisfied  the  equation 


H(r)  =  Va^+or^  -   a 

2 
where    a  -  M    -lo      This  means    the    shock  was   a*Bsumed   to 

be    a   hyperbolold    of  revolution,    with   a  radius    of 

curvature   equal    to    lo      The    a  was    so    chosen   that   the 

shock  approximates  a  sonic  disturbance  as  r  becomes 

largeo   This  condition  would  be  satisfied  by  every 

shock  created  by  the  motion  of  a  finite  blunt  body. 

Introduce  the  coordinates  x  and  y  defined  by: 
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X  =   z-H(r)  z   =  x+h(y) 

-  or  _      , 

7  -  r-^Q  r  -  y+r^ 

where   r      Is   a   constant    and  h(y)    =  H(Y+r    ). 
o  "  "^       o 

Since    the    flow   quantities    in   front    of   the    shock 
and  the    shock   shape    are    assumed,    the   flow  quantities 
immediately  behind  the    shpck  are    determined   as   functions 
of  y  by  the   Rankine-Hugoniot   shock  conditions   for  a 
polytropic   gas. 

eu2  („ 


(9+l)(l+h'^)+U^ 


Qp-1 

e-p 


u=U-      ^^-l^U 


e(i+h»^) 

V  =    (U-u)h' 

where  9  =  -^-^   ,    h'  stands  for  h'(y),  and  u,  v,  p,    and  p 
stand  for  u(0,y),  v(0,y),  p(0,y),  and  p(0,y)<,   Therefore 
the  expansion  of  u,  v,  p,  and  p  in  powers  of  y  can  be 
obtained  by  these  equations  from  the  expansion  of  h(y). 

The  flovj  quantities,  u,  v,  p,  p,  and  S  (entropy), 
satisfy  the  following  equations: 

(rpu)   +  (rpv)   =  0  (2) 

^^z  "^  ^"r  '^z^^   ^  °  ^^^ 

uS   +  vS   =  0  (5) 

z     r 

p  =  A(S)  pY  (6) 
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where  A(S.)  is  some  analytic  function  of  S,   (2)  is  the 
law  of  conservation  of  mass,  (3)  and  (I4.)  are  the  laws 
of  conservation  of  momentumj,  (5)  is  the  equation  which 
results  from  the  assumption  that  the  changes  are 
adiabatic,  and  (6)  is  the  equation  of  state  for  a 
polytropic  gaso   The  entropy  can  be  eliminated  between 
equations  (5)  and  (6)» 

These  equations  may  be  rewritten  in  the  x^j 
coordinate  system  in  the  following  form: 

S£  =  — 2_ . _______ — — (y) 

'^^  YP(l+ht-)-p(u-vh')^ 


h 


"Sx.  p(u  -  vh« 


(8) 


i£  „  h'  4h  +  pv  i^ 


Tx  p(u-vh«}  ^^^ 

1     p(u-vh« )4^  +  pv  4^  -  Ypv  4^ 

4e  =  1_ — _ii_^_E — iLJz  (10) 


The  only  term  which  could  possibly  not  be  analytic  is 

Yvp/(y+r  )  in  equation  (7)<.   If  this  were  singular 
o 

when  r   =  Oj,  it  would  be  impossible  to  compute  the  flow 
quantities  along  the  axis.   But,  presumably^  v  =  0  along 
the  axis;  hence  yvp/(y+r  )  should  be  analytic  and  give 
no  trouble  . 

In  terras  of  p(x,y.)  and  vCx^y);,  the  stream  function 
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^   is  given,  except  for  a  factor  of  2%,    by  the  equation 


^ 


=  ^|/(x,y)  =  I  U  (y+r^)^  -  (y+^o^  /  P(x',y)v(x',y)dx'. 


In  actual  computation  the  truncated  power  series  for 

J  p(x' ,y)  v(x%y  )dx'  is  obtained  by  multiplying  the 
o 
truncated  power  series  for  p  and  v  together  and 

integrating  term  by  termo 

Computation  of  the  Coefficients,   In  both  methods  of 

continuation  it  is  necessary  to  compute  the  coefficients 

of  the  truncated  power  series  of  the  unknowns  from  the 

differential  equations  and  the  initial  conditions «   In 

order  to  facilitate  the  discussion  of  the  computation 

of  the  coefficients  it  is  useful  to  introduce  the  following 

notation:   let  f.  .  be  the  coefficient  of  the  x  y^  term 

in  the  power  series  expansion  of  f(x,y);  that  is, 

i  1 
'  -^"^  where  the  sum  is  taken  over  all 


non-negative  i  and  j.   In  the  same  way,  let  h.  be  such 

J 

that  h(y)  =  ^^  h.y-^"., 

To  start  the  computation,  the  machine  reads  in  the 
four  constants  r  ,  k,  M,  and  y«   By  ^-^^  iteration  formula, 
the  machine  then  computes  h,  for  all  non-negative  J 
less  than  or  equal  to  k„   The  machine  then  uses  these 
coefficients  together  with  equations  (l)  to  compute  u   ., 
V      .,  p   .,  and  p   .  for  all  j  less  than  or  equal  to  k, 
and  places  them  in  appropriate  storage  locations.   This 
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ends  the  first  phase  of  the  computation.   The   computation 
of  the  coefficients  when  i  =  0  is  similar  to  the  computa- 
tion when  i  >  0  and  will  not  be  discussed  further. 

In  order  to  compute  the  coefficients  when  i  >  0, 
the  machine  contains  subroutines  which  can  form  the  sum, 
difference,  product,  and  quotient  of  two  power  series. 
To  be  explicit,  the  subroutine  is  instructed  on  entry  as 
to  which  coefficient  is  to  be  computed,  the  location  of 
the  coefficients  of  the  two  power  series  to  be  operated 
on,  and  the  location  where  the  coefficients  of  the  power 
series  of  the  sum,  difference,  etc.,  are  to  be  stored. 
There  are  also  subroutines  to  differentiate  and  integrate 
power  series.   The  subroutines  are  based  on  the  following 
simple  identities  (it  is  assumed  that  the  coefficients  of 
f(x,y)  and  g(x,y)  are  known): 
If  f(x,y)g(x,y)  =  h(x,y)  then 

I    J 

If  f(x,y)/g(x,y)  =  h(x,y)  then  h^  j  = 


^>^       ^  5^  ^i,j  ^I-1,J-J 


(11) 


I-l  J-1  J-1 


^0,0  ^'"^         i=0  j=0   ^*J   1-1, J-J    ^      1,J  o,J-j 
If  ^^^^f^  =   h(x,y)  then 
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If  f  (x,7)  =  ^^j^'^'^   and  J  >  1  then 

^I,J  =  1   ^I,J-1  •  ^^^^ 

In  (12)  If  either  I  =  0  or  J  =  0  the  appropriate  sum 

should  be  dropped. 

With  the  aid  of  the  subroutines  the  coefficients 

of  u,  V,  p,  and  p  for  i  >  0  are  easily  computed.   Assume 

that  u.  .,  V.  .,  p.  .,  and  p.  .  are  known  for  all 
i»J    ijJ    i»J       ^>J 

i+j  <  k,  i  <  I.   Prom  these  (4^)^  •,  ^J^^i  y    ®*^'  ^^® 
computed  for  i+j  <  k-1,  i  <  I  (differentiation  lowers 
the  order).   Every  power  series  on  the  right  of  equation 
(7)  is  known  for  i+j  <  k-1,  i  <  I,  hence  (-r^).  .  can  be 
computed  using  the  subroutines  for  i+j  <  k-1,  i  <  I. 
The  integration  subroutine  gives  p.  ,  for  i+j  <  k, 
i  =  I+l  (integration  raises  the  order).   In  a  similar 
manner  using  equations  (8)  to  (lO),  the  coefficients 
of  u,  V,  and  p  for  i+j  <  k,  1  =  I+l  are  computed.   Hence 
the  machine  essentially  does  an  induction  on  I  until  all 
of  the  coefficients  of  order  less  than  or  equal  to  k  are 
computed.   Actually  the  machine  stores  the  coefficients 
of  23  different  power  series  such  as  v  -j^  ,  p(u-vh'), 
etc.,  in  order  to  avoid  the  duplication  of  computations. 
Certain  restrictions  are  necessary  on  the  functions 
which  occur  in  the  differential  equations  if  actual 
computations  are  to  be  made  feasible.   Essentially,  the 
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functions  of  the  dependent  variables  must  be  easily 
expressible  as  power  series  in  the  independent  variables 
when  the  power  series  for  the  dependent  variables  are 
known  (i,e,  simple  recursion  formulas  Involving  the 
coefficients  of  the  power  series  of  the  dependent 
variables  must  give  the  coefficients  of  the  power  series 
of  the  functions )o   The  previous  section  has  just  shown 
that  functions  which  are  made  up  of  sums,  differences, 
products,  and  quotients  of  the  variables  and  their 
derivatives  are  of  this  type.   It  is  possible  for  the 
functions  to  contain  such  terms  as  sin(x)  or  even  V^  0 
but  sinCv),  for  example,  would  require  a  prohibitive 
amount  of  machine  computation  as  well  as  storage  space o 

Even  though  the  machine  used  was  the  IBM  70^,  the 
subroutines  are  essentially  those  used  by  the  Univac 
in  the  earlier  computations  done  by  R.  D.  Richtmyer. 
The  earlier  report  (6)  gives  a  more  detailed  account 
of  the  subroutines,  if  more  information  on  the 
subroutines  is  desiredo 

The  Entropy  Check..   Let  p,  p,  and  \j/  be  the  density, 
pressure,  and  the  stream  function  at  some  point  in  the 
flow.   The  streamline  through  this  point  presumably 
Intersects  the  shook.   Let  (0,y)  be  this  point  of 
intersection  in  the  x,y  coordinate  system  and  let 
p-,,  p^  be  the  density  and  pressure  on  the  downstream  side 
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of  the  intersection.   Since  the  stream  function  is  constant 
along  streamlines 

^  =  1  U(7+r^)2 

or  y  -  /TT  "  ^o  ' 

y   substituted   in  equations    (l)    gives   p-j^   and  p^.      Therefore, 

it    is   possible    to   determine   p.,    and  p,    from   the   knowledge 

of   ^. 

Since  entropy  is  constant  along  streamlines  it 

follows  from  (6)  that 

P_  =  (£_)Y 
Pi     Pi    " 
If  p,  p,  and  \!/  are  approximated,  it  is  possible  therefore 

to  find  an  approximation  to  (p/p-,  )  and  (p/p-,  )  o 

Comparing  these  last  two  approximations  will  give  an 

indication  of  the  accioracy  of  the  original  approximations. 

Input,  Output,  Etc    To  start  the  computation  r^,  k, 

Y  and  M  are  fed  into  the  m.achine  on  cards »   The  main 

code  then  computes  the  coefficients  of  the  power  series 

in  X  and  y  of  u,  v,  p,  p  and  pv  (used  to  compute  ^)    and 

punches  the  results  out  on  cards.   If  the  second  method 

of  continuation  is  being  used,  this  result  is  fed  back 

into  the  machine  with  the  constants  x  ,  y  ;  and  then, 

o'  "^o' 

the  coefficients  of  the  power  series  in  £,,    ^  of  u,  v, 
p,  p  and  pv  are  punched  out. 

With  the  first  method,  the  coefficients  of  the 
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power  series  in  x  and  y  are  fed  into  the  machine  viith 

the  output  code  plus  the  constants  x   and  y  .,   The 
^  -  0*^0 

machine  computes  the  power  series  in  E,   and  ^  and  reads 
in  a  control  cardo   With  the  second  method,  the  coeffici- 
ents of  the  power  series  in  E,   and  ^  are  fed  into  the 
machine  with  the  output  code  and  a  control  cardo 

This  control  card  for  both  methods  contains  [skj, 
^l  ,  6C,  ^,    ^-1,  &K>,    and  no   These  last  six  constants  are 
to  instruct  the  nachine  to  print  out  on  the  tabulator 
the  flow  quantities  at  the  points  (4, +154,  Z-,'^2^K) 
where  i  =0  to  m  and  j  =  0  to  n.   First  the  machine 
prints  all  the  flow  quantities  at  the  points  (^-j  +  lS^,^-) 
i  =  0  to  m,  and  then  proceeds  to  search  for  the  point 
on  the  surface  of  the  body  with  ^  =  ^, »   In  order  to 
do  this,  it  is  necessary  to  solve  the  equation  \1/ (4,  ^-.  )=0, 
This  is  done  by  Newton's  method  for  finding  the  roots  of 
an  equation^   The  flow  quantities  are  then  printed  out 
for  this  point  on  the  body.   Next  the  machine  prints  out 
the  flow  quantities  for  the  points  i£,^+16^gil,-.+6^) 
i  =  0  to  mj  etco 

The  following  flow  quantities  are  printed  out: 
^}    V,  p,  p,  \1/,  (p/p-j^)-,  (p/p-|_),  and  (u  +v  )o   (u  +v  )  is 
printed  out  in  order  to  have  a  simple  indication  as  to 
whether  the  flow  is  stibsonic  or  supersonic »   Let   q  be 
the  Bernoulli  number  and  c  be  the  critical  speed, 
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These  are  related  by  the  equation: 

4:.2  _  Y-i  a2 

(For  Y  =  l.k-   and  M  =  12,  c'*^^  =  3^.766666.) 

p   p     ''*p 
If  (u  +v  )  >  c"   the  flow  is  supersonic. 

If  (u  +v  )  <  c"   the  flow  is  subsonic. 
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Appendix  Bo   The  Arithmetic  Used  in  the  Testo 

V/hen  the  power  series  method  for  solving  the 
detached  shock  problem  was  first  coded  for  the  Unlvac, 
a  straight  forward  single-precision  floating-decimal 
arithmetic  was  usedo   This  arithmetic  proved  unsuccessful. 
Some  of  the  high  order  coefficients  were  so  Inaccurate 
that  it  was  necessary  to  discard  these  coefficients 
before  the  power  series  could  be  evaluated,,   In  order 
to  do  this  it  was  necessary  to  distinguish  between  the 
accTorate  and  the  inaccurate  coefficients o   This  was 
accomplished  by  receding  the  detached  shock  problem 
with  single-precision  floating-decimal  subroutines  with 
significance  monltoringo   Each  numerical  quantity  was 
represented  by  not  only  a  fraction  and  an  exponent  but 
also  a  significance  index o   This  significance  index  was 
an  integer  which  indicated  roughly  how  many  digits  in 
the  fraction  were  accurate.   This  monitoring  scheme 
proved  very  successful,   (See  [6]  for  a  description  of 
these  subroutines,,) 

When  the  problem  was  receded  for  the  IBM  70)4-  it 
was  decided  to  use  double-precision  and  to  retain 
accuracy  monitoring  but  to  use  an  alternative  procedure 
suggested  by  No  Co  Metropolis »   Ordinarily  the  fractional 
part  of  a  floating-point  number  on  the  7O-!-  (a  binary 
machine)  has  a  leading  oneo   Instead  of  keeping  a 
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significance  index,  enough  leading  zeros  are  supplied 
to  the  fraction  so  that  all  the  remaining  bits  are 
significant.   Because  leading  zeros  are  allowed  in  the 
fraction,  this  arithmetic  is  called  unnormalized.   This 
monitoring  scheme  has  the  advantage  that  no  bits  are 
wasted  in  storing  the  significance  index. 

Accuracy  monitoring  for  this  unnormalized  floating- 
point arithmetic  Is  done  as  follows:   Ordinarly,  when 
two  floating-point  numbers  are  added  or  subtracted,  the 
last  step  in  the  operation  consists  in  normalizing  the 
answer  (that  is,  the  fraction  of  the  answer  is  shifted 
left  and  the  exponent  is  reduced  until  there  is  a  leading 
one  in  the  fraction).   In  the  present  arithmetic  this 
last  step  is  suppressed.   In  multiplication  and  division 
the  niOTiber  of  leading  -zeros  in  the  fraction  of  the  answer 
is  exactly  the  same  as  the  number  of  leading  zeros  in 
the  least  significant  of  the  two  n"umbers  entering  into 
the  operation. 

To  check  the  monitoring,  the  arithmetic  was  used 
to  compute  sin(x)  by  its  power  series  for  large  values 
of  X.   The  results  are  listed  in  Table  l^.,  page  67* 
Column  A  lists  the  number  of  significant  digits  (i.e. 
not  counting  the  leading  zeros)  in  the  fraction  out  of 
the  possible  Sk-'      Column  B  lists  the  nxjmber  of  accurate 
bits  in  the  fraction.   Column  C  lists  the  number  of 
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accurate  bits  in  the  fraction  when  normalized  arithmetic 
is  used»   Considering  that  the  monitoring  scheme  does 
not  take  into  account  the  accumulation  of  error  due  to 
round-off,  the  results  are  surprisingly  good. 

Details  of  this  arithmetic  and  the  code  itself 
may  be  obtained  by  writing  to  the  SHARE  Librarian, 
590  Madison  Avenue,  New  York  22,  New  York,  and  asking 
for  NU-UDPlo 

In  all  of  the  computations  so  far  the  significance 
monitoring  was  Indispensable o   Even  though  the  coefficients 
computed  on  the  7OI4-  usually  have  a  slnglf  icance  of  30 
bits  or  more,  no  reliance  could  be  placed  on  the 
coefficients  of  large  order  without  the  monltoringo 
Now  that  more  work  has  been  done  on  the  shock  problem. 
It  is  no  do\±>  t  possible  to  eliminate  the  monltoringo 
This  would  mean  a  great  saving  in  time.   Since ^  in  all 
of  the  computations,  the  significance  of  the  coefficients 
varies  inversely  x^rith  the  order,  it  would  seem  reasonable 
to  rely  on  this  fact  in  future  computationSe   k   can  be 
varied  and  the  accuracy  of  the  computation  can  be 
checked  (by  means  of  Bernoulli ''s  law  and  the  entropy 
check)  until  a  reasonable  k  is  foundo   If  k  is  too 
large  the  approximation  is  corrupted  by  round-off  error; 
if  k  is  too  small  the  approximation  is  corrupted  by 
truncation  error o 
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Table   1 
The    Location   of    the    Body 
Results   Using  the   Second  Method   of   Continuation 


r 

z 

u 

V 

P 

P 

0 

.10153518 

0 

0 

6.237223 

185.8692 

.05 

.10328060 

.039276 

.5623781 

6.213i;93 

181;. 8800 

.1 

.  10852381; 

.155788 

1.113109 

6.114-3220 

181.9593 

.15 

.11728574 

.3I1-5687 

lo6[;1201 

6.029069 

177  0  214-31; 

.2 

.1296009 

.602890 

2„1368V7 

i5. 87629 

170. 9Mi? 

.25 

.1[;55182 

.919511 

2.59213 

5.68708 

163.329 

.3 

.I6!50y96l; 

1.28661 

3.000793 

5.1;7065 

151;.691; 

.35 

.1881^216 

l,69i;3 

3.35889 

S.2.1ZS 

11;5.35 

.1+ 

.21613291; 

2.132750 

3.66)|J|11; 

ij.o9790l;8 

135.5891 

.M-5 

.2ii666008 

2.592515 

3,917230 

I;c7l6305 

125 06789 

.5 

.2817980 

3.06i;71; 

1;<.11878 

i;.)|)l98l 

115.81;9 

SS 

.321123 

3.5i;00 

1;.2722 

[;.1858 

106.33 

.6 

06k8i; 

3.988 

1;.388 

3.9I;9 

97.87 

.7  .1+6570233  i;.91;1076  );.l;762l;8 

.75  .S2'}>'5S'iSZ  5.382131  l;.I;7i;31;9 

.8  .5860619  5.80537  l;.Wl-83 

.85  .651l066  6.2057  4.3932 

.9  .72829  6.528  1;.31;9 


3<.1;261057  80. 34169 

3.194297  72.83585 

2.97533  65c 9419 

2.7737  59.739 

2.667  55.88 
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Table  2 
Results  Using  the  First  Method  of  Continuation  on 
Elementary  Functions  with  x-,  =  1.2 

Description 

The  results  are  tabulated  (on  the  next  page) 
when  the  first  method  of  continuation  is  used  to  compute 
the  value  of  the  functions  l/(l+x),  l/(l+x)^,  l/(l+x^), 
and  l/(l+x-^)  at  X,  =1.2  when  only  the  first  k  terms 
of  the  power  series  about  the  origin  of  these  functions 
are  used.   In  accordance  with  the  discussion  on  page  1^2, 
X   and  s  were  chosen  in  the  following  way: 

X  =  .I4.1392785  and  s  =  .235  when  computing 


l/(l+x)  and  l/(l+x)^. 


X      =    .51720723    and   s   =    .202  when    computing   l/(l+x    ). 
X      =    ,62752519    and   s    =    .ll|7  when   computing   l/(l-i-x-^/. 

The    columns   headed    [ sk ]    are    actually    [sk+1/2]. 

The    error   listed   is   the    approximation  minus   the   true 

value.      8.    -3    should  be   read   8   times   10      . 
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Results  Using  the  First  Method  of  Continuation 
on  Elementary  Functions  with  x,  =  lo2 

The  description  of  the  table  is  on  the  previous  page, 


l/(l+x) 

l/(l+x)^ 

l/( 

1^x2; 

1 

l/(i+x^: 

) 

k 

[sk] 

error 

[sk] 

error 

[sk] 

error 

[sk] 

error 

0 

0 

5c 

-1 

0 

9o 

-1 

0 

6» 

-1 

0 

6. 

-1 

5 

1 

-2. 

-1 

1 

-1. 

0 

1 

7o 

-2 

1 

-3. 

-1 

10 

3 

1. 

-2 

3 

l|o 

0 

2 

-1. 

-1 

1 

-3e 

-2 

15 

i^. 

1, 

-2 

h 

-2  = 

0 

3 

5e 

-2 

2 

-2c 

-1 

20 

5 

-1, 

-2 

5 

7» 

-1 

1+ 

3o 

-3 

3 

7o 

-2 

2.S 

6 

7o 

-3 

6 

-2o 

-1 

5 

-2. 

-2 

1+ 

Ic 

-1 

30 

8 

i|» 

-3 

8 

8o 

-1 

6 

Ic 

-2 

1+ 

7o 

-2 

35 

9 

-2. 

-3 

9 

-3o 

-1 

7 

-2  c 

-3 

5 

-3. 

-2 

lj.0 

10 

8„ 

-i| 

10 

1. 

-1 

8 

-5» 

-3 

6 

-3. 

-2 

hS 

11 

-^. 

-1^ 

11 

-5o 

-2 

9 

l+o 

-3 

7 

8c 

-3 

50 

13 

-8. 

-5 

13 

5. 

-2 

10 

lo 

-1+ 

7 

2, 

-2 

SS 

\h 

6, 

-5 

11; 

-2c 

-2 

11 

-2o 

-3 

8 

Ic 

-2 

60 

15 

-3c 

-5 

15 

1^-c 

-3 

12 

8c 

-1+ 

9 

-6o 

-3 

65 

17 

-2c 

-5 

17 

-lo 

-2 

13 

3c 

-I; 

10 

-5o 

-3 

70 

18 

8. 

-6 

18 

5. 

-3 

11+ 

-[+c 

-i| 

10 

-2  c 

-3 

75 

19 

-1|. 

-6 

19 

-2c 

-3 

15 

1, 

-1+ 

11 

3. 

-3 

80 

20 

2c 

-6 

20 

7. 

-i; 

16 

6c 

-5 

12 

.9o 

-1+ 

85 

22 

l+o 

-7 

22 

-7o 

»[|. 

17 

-7. 

-5 

12 

2c 

-1+ 

90 

23 

-3o 

-7 

23 

2c 

-1^- 

18 

2. 

-5 

13 

-1. 

-3 

95 

2i^. 

2o 

-7 

21+ 

-3. 

-5 

19 

1, 

"5 

111. 

1|. 

-5 

100 

26 

9, 

-8 

26 

Ic 

-1; 
k\  - 

20 

-2, 

-5 

15 

1+. 

-1+ 

Table  3 

Results  Using  the  First  Method  of  Continuation  on 

Elementary  Functions  with  x,  =  .8 

Description 

The  results  are  tabulated  (on  the  next  page) 
when  the  first  method  of  continuation  is  used  to  compute 
the  value  of  the  functions  l/(l+x),  V   (1+x  ),  and 
l/{l+x-^J  at  x,  =  .8  when  only  the  first  k  terms  of 
the  power  series  about  the  origin  of  these  functions 
are  usedo  'In  accordance  with  the  discussion  on  page  \\.2, 
X   and  s  were  chosen  in  the  following  way: 

X  =  .2I4.829866  and  s  =  .14.38  when  computing  l/(l+x). 
X   =  <,296L|.8953  and  s  =  «l\2S   when  computing  l/(l+x  ). 
X   =  .314-014.3853  and  s  =  .1|18  when  computing  l/(l+x-^). 

The  columns  headed  [sk]  are  actually  [sk+l/2]. 
Under  each  function  two  errors  are  listed.   The  first 
error  coliAmn  is  the  approximation  minus  the  true  value 
when  the  first  method  of  continuation  is  used.   The 
second  error  column  is  the  approximation  minus  the  true 
value  when  no  additional  truncation  is  used  (i.e.  s  =  l) 
8.  -3  should  be  read  8  times  10  -^. 
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Results  Using  the  First  Method  of  Continuation 

on  Elementary  Functions  with  x,  =  08 

The  description  of  the  table  is  on  the  previous  page. 

l/(l+x)          l/(l+x^)  l/(l+x^) 

k    [sk]  error  error   [sk]  error  error  [sk]  error  error 

0    0  1^,  -1  ko   -1  0  i+«  -1  k-   -1  0  3o  -1  3.  -1 

5     2   3.  -2  -1.  -1  2  -60  -2   2.  "1  2  -80  -2  -2o  -1 

10     k     lo  -2  5o  -2  k     2o  -2  -ij..  -2  [4.  -6c  -3  -5.  -2 

15     7  -1«  -3  -2o  ~2  6  -3.  "3  -2,  -2  6   3^  -3  -1.  "2 

20     9  -8.  -5  5o  -3  8  -8«  -5  5.  -3  8  -l.  -3  60  -3 

25    11  -[;,  -5  -2.  -3  11  1.  -i;  2,  -3  10  6«  -1|  2o  -3 

30    13  -5.  -6  6.  -k  13  -60,  -6  -5o  "Ii-  13  5.  -6  U-"  -k- 

3^    15  -1.  -6  -2.  -ij.  15  -6»  -6  -2„  -1|  15  "1.  -6  -2,  -14. 

l|.0    18  2.  -7  6„  -5  17  2.  "6  5o  -5  17  -60  -6  -6.  -5 

1^5    20  l„  -8  -2.  -5  19  -3»  -7  2o  -5  19  1.  -6  -1.  -5 

50    22  5„  -9  6„  ~6  21  6.  -8  -60  -6  21  -3o  -7  8„  -6 

55    2I(.  60-10  -2.  -6  23  -7o  -9  -2o  -6  23  1-  -7  2.  -6 

60    26  2„-10  7.  -7  25  -1»  -9  60  -7  25  -2c  -8  5o  -7 

65    28  3.-11  -2.  -7  28  1^.-10  2o  -7  27  2.  -9  -3«  -7 

70    31  -2.-12   7.  -8  30  6.-12  -60  -8  29  -5o-10  -7o  -8 

75    33  -60-13  -2„  -8  32  -3.-11  -3o  -8  31  "2.-10  -2.  -8 

80    35  -8. -Ill  8„  -9  3k     7.-12  7<.  -9  33  lo-lO  9.  -9 

85    37  -2.-1!+  -3.  -9  36  -1.-12   3^  -9  36  2.-11  2o  -9 

90    39  -Ivo-15  80-10  38  2o-13  -7.-10  38  5o-13  6.-10 

95    1|2  -6.-16  -3.-10  i;0  -1.-11+  -3.-10  1+0  -1.-12  -3.-10 

100    I44  -7. -16  9.-11  [^2  -9.-15   8.-11  1+2  3o-13  -8.-11 
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Table  I4. 

Unnormazlied  Arithmetic  as  an  Accuracy  Monitor 
Computation  of  sm  (j-  +  ^nj 

Column  A  is  the  number  of  bits  in  the  sine  using 

unnormalized  arithmetico 

Coltimn  B  is  the  number  of  accurate  bits  in  the  sine 

using  unnormalized  arithmetic. 

Column  C  is  the  number  of  accurate  bits  in  the  sine 

using  normalized  arithmetic* 

■      n         A         B         C 


0 

.51+ 

53 

53 

1 

52 

50 

50 

2 

50 

1+8 

51+ 

3 

k8 

1^7 

1+7 

1^ 

k.6 

hh 

hk 

5 

hk 

k2 

k6 

6 

k2 

1^0 

39 

7 

ko 

39 

1+0 

8 

37 

3B 

31+ 

9 

35 

35 

33 

10 

33 

32 

3k 

11 

31 

30 

29 

12 

29 

28 

29 

13 

26 

26 

2S 

1^    ■ 

21| 

21 

■21 

15 

22 

19 

20 

16 

20 

20 

•     17 

17 

18 

18 

17 

18 

15 

13 

17 

19 

13 

10 

11 

20 

11 

■  10 

10 

21 

9 

9 

9 

22 

6 

1+ 

1+ 

23 

k 

3 

k 

21+ 

0 

0 

h 
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